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ABSTRACT 


This dissertation focuses on the signal processing problems associated with the 
detection of hazardous windshears using airborne Doppler radar when weak weather 
returns are in the presence of strong clutter returns. In light of the frequent inade- 
quacy of spectral-processing oriented clutter suppression methods, we model a clutter 
signal as multiple sinusoids plus Gaussian noise, and propose adaptive filtering ap- 
proaches that better capture the temporal characteristics of the signal process. This 
idea leads to two research topics in signal processing: (1) signal modeling and pa- 
rameter estimation, and (2) adaptive filtering in this particular signal environment. 
A high-resolution, low SNR threshold ML frequency estimation and signal modeling 
algorithm is devised and proves capable of delineating both the spectral and tem- 
poral nature of the clutter return. Furthermore, the LMS-based adaptive filter’s 
performance for the proposed signal model is investigated, and promising simulation 
results have testified to its potential for clutter rejection leading to more accurate 
estimation of windspeed thus obtaining a better assessment of the windshear hazard. 
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CHAPTER 1 


INTRODUCTION 
1 . 1 Modeling Discrete-Time Signal 

Though conventional time series (signal) analysis is heavily dependent on the 
twin assumptions of linearity and stationarity, since the late 1960s, parametric mod- 
eling of nonstationary signals has received a great deal of attention [1]. However, the 
definitions of signal nonstationarity are still diversified among engineering researchers 
and mathematical statisticians. As opposed to statisticians’ interest in the tempo- 
ral characteristics of an observed data record, signal processors are generally more 
absorbed in analyzing its spectral content [2, 3]. In Section 1.1.1, we give a brief intro- 
duction of some commonly used time-domain approaches for modeling nonstationary 
signals. 

1.1.1 Time Series Analysis Approach 

It is generally recognized by statisticians that nonstationary processes may arise 
in several ways [1, 4]. 

The first is the “trend plus stationary residual” model 

x(n) = n(n) + a(n) 

where n{n) is a deterministic function, and a(n) is a zero-mean stationary processes. 
An extension of this model, i.e. , the “trend-seasonal-irregular” model, 

x(n) = T[n) -f S(n) + R(n) 

with T(n ) as the trend, 5'(n) the seasonal term, and R(n) an irregular component, is 
often encountered in econometric time series analysis [5]. 

The second is the the autoregressive integrated moving average (ARIMA) process 

MB){ 1 - B) d x(n) = 9 q (B)e{n) 



where B is a backshift operator such that B k x(n) = x(n — k ), the stationary AR 
operator <f> p (B) = 1 — a\B — ••• — a p B p and the invertible MA operator 0 q (B) = 

1 — b\B b q B q share no common factor, and e(n) is a white noise process of zero 

mean and constant variance. This model postulates that differencing the x(n) process 
d times will result in a stationary ARMA process. Furthermore, one should notice 
that the nonstationary characteristics of this model may be seen by interpreting x(n) 
as the result of passing the white noise process e(n) through an “unstable” HR filter 
whose transfer function — — 777 has d overlapping poles on the unit circle. 

- z ~ x Y 

The third is the ARMA process with freely varying time-dependent parameters 
which is a more general class of nonstationary processes in which the time-dependent 
nature of the processes can be delineated in various ways [1, 6]. 

Many common signals analyzed in practice are indeed not stationary, and their 
time-domain characteristics may be captured by the models proposed above. How- 
ever, in the development of various spectrum estimation algorithms, as addressed by 
Marple [ 2 ] and Kay [ 3 ], short data segments from the longer data record may be con- 
sidered to be locally stationary , availing those estimation methods to many real world 
applications. Assuming the data record x(n) is from a stationary process, the formal 
definition of the spectrum shows that it is a function strictly of the second-order 
statistics [2]. The second-order statistics are also assumed to remain unchanged, or 
stationary , over time. Thus, the spectrum is not a complete statistical picture of a 
random process that may have other information in third- and higher-order statistics. 
Section 1 . 1.2 deals with another modeling approach particularly suitable for a nar- 
rowband signal process, where the signal’s features over both the time-domain and 
the frequency- domain are simultaneously retained. 

1.1.2 Sinusoidal Parameter Estimation Approach 
The estimation of the frequencies of sinusoidal components embedded in white 
noise is a problem that arises in many fields such as communications, radar, and sonar 
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applications, where enhancement (suppression) of narrowband signals (interferences) 

is among the signal processors’ major concerns [3, 7]. In this model, given as 

v 

x(n ) = ^2 Ai exp(j2irfin + <f>i) + w(n) 

»=i 

where w(n ) is complex white Gaussian noise of zero mean and constant variance 
[8], the signal amplitudes are usually assumed to be constants, whereas the phases 
are either constants or random variables dependent on the frequency estimation ap- 
proaches: the approximate maximum likelihood techniques regard the sinusoidal pro- 
cess as a deterministic signal with unknown frequencies, while the eigenanalysis ap- 
proaches employ a WSS random process model so that the frequencies appear as 
unknown parameters in an autocorrelation matrix. It should be pointed out that 
when the sinusoids are of constant phase, this model belongs to the first category of 
nonstationary processes discussed in the previous section. 

As will be addressed in next two chapters, this dissertation focuses on the maxi- 
mum likelihood frequency estimation problems with the idea of modeling the observed 
data as a deterministic sinusoidal signal with unknown parameters plus noise, though 
the data itself may be a partial realization of a nonstationary process [3, 9]. 

1.2 Motivation for Adaptive Processing Scheme 
One of the most challenging problems facing modern signal processors is the de- 
tection and estimation of weather information using airborne Doppler radar [10, 11]. 
This application is associated with the development of a new generation airborne 
Doppler weather radar for detecting windshear which is particularly hazardous to 
aircraft operating at low altitudes [11, 12, 13]. For an aircraft in a landing mode, 
with an antenna scanning the airspace along the intended flight path, the radar re- 
turn represents a sample function of a nonstationary process which in the presence 
of weather can consist of two slowly time-varying components. A signal component 
is due to the weather return and a noise component is due to the ground clutter re- 
turn. Furthermore, in situations where the signal-to-clutter ratio (SCR) is extremely 
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low, it becomes a very difficult task to suppress the clutter component so that the 
radar return can be processed to determine windfield characteristics within the search 
volume. 

The general approach in reducing the effects of ground clutter on windspeed 
estimates made from the Doppler radar return sequences is to use some form of 
frequency selective filtering [10]. Fixed-notch clutter rejection filters are commonly 
used to separate the two processes on the basis of spectral content. Moreover, in most 
spectral analysis oriented approaches, two inherent problems are unavoidable. First, 
good spectrum estimation can only be attained when the signal-to-noise ratio (SNR) 
is above a certain threshold; Secondly, the signal’s temporal information (e.g., phase) 
is often discarded. Accordingly, approaches of this sort will fail to enhance the weak 
weather return in the presence of strong clutter returns. 

An adaptive filter, on the other hand, is an active system that will make appro- 
priate parametric adjustments to the unknown signal environment or its time- varying 
temporal statistics, by following some type of optimization criterion. For the applica- 
tion of clutter rejection and weather information extraction associated with airborne 
Doppler radar, some promising results obtained by LMS-based adaptive filtering have 
been reported by Lai and Baxa [14]. This dissertation serves as an in-depth study of 
their work. Major contributions are summarized in the following section. 

1.3 Contribution and Organization of This Dissertation 

The contributions of this dissertation include three aspects. First, as noted 
in Chapter 2, the major difficulty in ML frequency estimation is that the likelihood 
function is highly nonlinear with respect to the frequency parameters. Moreover, given 
a short data record under low SNR conditions, the presence of spectrally close signals 
makes correct detection a very tough task. Chapter 2 formulates this ML frequency 
estimation problem, briefly surveys its solutions proposed by the signal processing 
community, and provides the basic ideas in Wiener filtering and the adaptive LMS 
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algorithm. All the discussions therein serve as a theoretical basis for the chapters to 
follow. 

In Chapter 3, we devise a high-resolution frequency estimation and signal detec- 
tion scheme utilizing the expectation-maximization (EM) algorithm. This method 
efficiently solves the problem of multi-dimensional optimization, and through the in- 
clusion of MDL criterion or Fisher’s T 0 statistic, attains high probability of correct 
detection in low SNR situations. 

Secondly, convergence responses of the LMS-based adaptive noise canceler to 
sinusoidal signal processes are investigated in Chapter 4, through extensive simula- 
tion. Though the approach taken is somewhat more experimental than analytical, 
our simulation results provide valuable insights regarding different adaptive filter 
configurations and design considerations, under a very realistic signal environment 
encountered in many applications. 

Thirdly in Chapter 5 the adaptive filter is investigated through simulation as 
a process decorrelator to suppress the clutter return and thus enhance the weather 
return. Results from Chapter 4 have been used as a design guideline for these adaptive 
filtering schemes, whereas the algorithms in Chapter 3 are employed to model the 
Doppler weather radar clutter return process. As our simulation results show, even 
in a low SCR signal scenario, an adaptive signal processing scheme can significantly 
outperform traditional fixed filtering approaches. In light of all these analysis and 
simulation results, a concise conclusion and some suggestions for future work are 
provided in Chapter 6. 
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CHAPTER 2 


ML FREQUENCY ESTIMATION AND ADAPTIVE FILTER 

In this chapter, we give a brief overview of two important research topics in signal 
processing: spectral estimation and adaptive filtering. All the discussions here will 
serve as the theoretical basis for Chapters 3 and 4 and their applications in Chapter 
5. Traditionally, power spectral density estimation has very much relied upon the 
Fourier transform, the theories of random processes and filter theory. Specifically, 
classical estimation methods like periodogram and the Blackman- Tukey spectral es- 
timator are Fourier-transform-based, whereas “modern” spectral estimation, which is 
a parametric modeling approach, mainly depends on random process (e.g., ARMA) 
and filter theories [2, 3]. Our attention will be directed to the problem of parameter 
estimation of sinusoids corrupted by noise, a research issue still challenging many 
signal processors, especially when a short data record is available and high frequency 
resolution is required [3, 15, 16]. An adaptive filter is an active system that will 
make appropriate parametric adjustments to the unknown signal environment or its 
time-varying statistics, by following some type of optimization criterion [7, 17]. Ac- 
cording to Haykin [7], to derive the recursive algorithms for the operation of adaptive 
filters, one can identify three distinct methods: (1) an approach based on Wiener 
filter theory, (2) an approach based on Kalman filter theory, and (3) the method of 
least squares. We will focus on the first approach, Wiener filtering. 

2.1 Basic Assumptions in Signal Model 
To begin the discussion of maximum likelihood estimation (MLE) of frequency 
parameters, specifically assume that the received data vector y = [y (0) j/( 1 ) • • • y(N— 

1)] T consists of p complex sinusoids in complex white Gaussian noise (CWGN) 

v 

y( n ) = ]C A ' ex P(j’27r/.n) + w(n) 

t=l 


( 2 . 1 ) 



where A, = \Ai\e j4> ' is the complex amplitude of the ith sinusoidal component, and 
w(n) is complex white Gaussian noise with zero mean and variance a 2 . The sinu- 
soidal parameters {|Ai|, l^l, • • • , |A P |, <j>i, <i> 2 , ■ • * , /21 •••>/?}> which consists of 

amplitudes (0 < |A,| < 00 ), phases (0 < <j>i < 27 r), and frequencies (0 < /, < 1), are 
assumed to be constant but unknown and are to be estimated. The idea here is to 
model the observed data as a deterministic sinusoidal signal with unknown param- 
eters plus white noise, though the data itself may be a snapshot of a nonstationary 
process [3, 9]. It is also important to remember that the discussion to follow assumes 
that the number of sinusoids, p, is known. Methods to estimate p are discussed in 
Section 2.3. 


2.2 Problem Formulation and Implementation of MLE 


2.2.1 Single Sinusoid 

To help formulate the ML estimation problem, let us first look at the case of 
a single sinusoid. Let ej = [1 e j2irfl e j2 * fl2 ••• e i 2 *M N ~ l )] T . It is well known in the 
literature that the MLE of the frequency and amplitude of a single complex sinu- 
soid in complex white Gaussian noise is found by minimizing the scoring function , 
S(Ai,fi) = ||y - Aiei|| 2 = (y - ^iei) H (y - Aiei), where H denotes the Hermitian 
transposition and || • || the Euclidean norm. This is equivalent to choosing the fre- 
quency at which the periodogram attains its maximum [3]. The MLE of { |v4i | , <^ 1 , / 1 } 
is given as 


h 


1 

argm /fiv 



1 

arg max — 

h N 


N—l 

Y, y{n)exp(-j2Trf l n ) 


n=0 


2 


jV-1 


I Ail = 

<^1 = arctan 


i H y(n)exp(-;'27r/ 1 n) 

n=0 


( 2 . 2 ) 


Im (E^o 1 y(») exp(-j27r/ 1 n) 

.Re (EnJo 1 J/(n)exp(-j27r/ 1 n) 

Note that the frequency estimate is found as the result of a one-dimensional search. 
To achieve this goal, a coarse FFT followed by a Newton-Raphson search can be 
implemented. Please refer to the Appendix for details. 
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2.2.2 Multiple Sinusoids 

In vector form, the observed data y = [?/(0)y(l) • • • y(N— 1)] T can be represented 


y = Ai e, 1 + w = EA + w (2.3) 

t=i 

where e, = [1 e j2irf ' e j2 * fi2 ••• e^ 2ir ^ N ~^] T , E = [ej e 2 ---e p ], A = [A\ A 2 • • • A P ] T , 
and w = [tu(0) u>(2) • • • w(N — 1)] T . For notational simplicity we have concealed the 
dependence of matrix E on the frequency parameter vector f = [fi / 2 • ■ • f p \. 

Based on the white Gaussian assumption, it has been shown that the MLE of 
A and f, denoted as A and f, can be found by solving the following nonlinear least 
square problem [3]: 


{A,f} 

S(A,f) 


Fixing f and thus assuming 


= arg min 5(A, f) 

= ||y-EAf = (y-EA)"(y-EA) 


N - 1 


n=0 


y{n) - ^2 A. exp(;'27r/,n) 


1=1 


(2.4) 


(2.5) 


E to be a known matrix, S( A, f ) is minimized over A by 


A = (E^E) - ^^. 


(2.6) 


It is well known that if f is replaced by its MLE f, then A will be the MLE of A [3]. 
Substituting (2.6) into (2.5) yields the scoring function 

S( A,f) = y H (y-EA) 

= y"y - y"E(E H E)- 1 E"y. (2-7) 

Therefore, the MLE f of the frequencies can be found as 

f = arg max L(y ; f ) 

where L(y;f) = y H E(E // E)" 1 E i/ y. (2.8) 


L( y; f ) is a highly nonlinear function of the unknown frequencies and therein lies the 
central problem of ML estimation. We refer to L( y; f) as the log-likelihood of y given 
f since it is In f y (y, f ) with terms not associated with f discarded. 
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2.2.3 Implementation of MLE 


Modern methods for ML frequency estimation generally fall within two classes: 
linear prediction based methods and eigenstructure-based methods [3, 16]. In the 
' 1980s, the modified forward-backward linear prediction (MFBLP) approach developed 
by Tufts and Kumaresan [15] achieves great success in that it outperforms other 
eigendecomposition-based methods in terms of its good frequency resolution at low 
SNR conditions. 

However, since Ziskind and Wax [9] proposed the alternating projection (AP) 
algorithm in the late 1980s, several AP-oriented approaches that show even better 
performances than that of the MFBLP method have been developed by signal proces- 
sors [16, 18, 19]. In Chapter 3, we will devise a new class of high- resolution frequency 
estimation algorithm that not only possesses some nice computational and statis- 
tical features, but also demonstrates performances favorably compared with those 
advanced methods. 


2.3 Order Selection 

In this section, we consider two information theoretic criteria [20] found favorable 
in both time series analysis and signal processing for model order selection, and give 
their modified versions for deciding the number of sinusoids buried in complex white 
Gaussian noise. Let /(y ; 0) denote the PDF of data y given the true parameter vector 
0 of p components. Assuming we have an MLE 0 of 0 for each model of order k (i.e., 
the number of independently adjusted parameters in the model) in question, Akaike’s 
information criterion (AIC) [21] is defined as 

AIC(k) = -2 In /( y; 0) + 2Jb, (2.9) 

and according to Wax and Kailath [20], the minimum description length (MDL) 
principle originally proposed by Rissanen [22] is given as 

MDL(k) = - In /(y; 0) + ]-k In N. (2.10) 
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To tailor AIC for modeling the signal described by (2.3), Kay [3] suggests that we 
can choose the number of sinusoids by minimizing 

AIC (i) = 2N In Si{ A J)+6i (2.11) 

where i is the number of sinusoids, and k in (2.9) has been replaced by 3i to account 
for the three parameters (amplitude, phase, and frequency) associated with each 
sinusoid. Before using MDL, we need to remember that Rissanen’s basic idea behind 
this order selection principle is to find the shortest code length (number of bits) to 
encode an observed data set y, assuming the parameter 0 is a vector of k real-valued 
components [22]. As expressed in (2.5), the parameters to be estimated are associated 
with the complex-valued multiple sinusoidal signal vector EA. Therefore, replacing k 
in (2.10) by 2k, and following the same rationale used by Kay to account for the three 
parameters associated with each sinusoid, the MDL criterion can be presented as 

MDL(i) = N In 5,( A, f ) + 3i In N (2.12) 

where i is the number of sinusoids. Notice that the first terms in (2.11) and (2.12) 
involve the negative of the log-likelihood of an ML estimate, and the second terms 
represent the “cost” for the selected model order (correspondingly the number of 
parameters). However, as observed in (2.12), the MDL criterion will “penalize” more 
for redundancy in modeling than the AIC. Performance comparison of these two 
model-order estimators will be given in Section 3.4 through simulation. 

2.4 Wiener Filtering and Widrow-Hoff LMS Algorithm 
Adaptive filters are mainly derived from linear optimum filter theory for wide- 
sense stationary stochastic processes [7, 17]. Regarding the filter specification, two 
choices have to be made clear. First, the choice of a finite impulse response (FIR) or 
an infinite impulse response (HR) for the filter is dictated by practical considerations. 
Secondly, the type of statistical criterion used for the optimization is often influenced 
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by mathematical tractability [7]. Our attention will be confined to the more well- 
developed FIR filter theory and the least-mean-square (LMS) algorithm, a simplified 
criterion derived from the method of steepest descent. 

2.4.1 Wiener Filter 

This section begins with the discussion of a class of optimum linear discrete-time 
filters known collectively as Wiener filters. Consider the block diagram of Figure 2.1 
for this specific filtering problem. The filter input consists of a time series y( 0), y(l), 
y(2) • • •, and the filter is itself characterized by the impulse response w 0 , w t , w 2 
At some time instant n, the filter produces an output denoted by u(ra). This output is 
used to provide an estimate of a desired response d(n). With the filter input and the 
desired response representing single realizations of respective stochastic processes, the 
estimation is accompanied by an error with statistical characteristics of its own. In 
particular, the estimation error , denoted by e(n), is defined as the difference between 
the desired response d(n) and the filter output u(n). The basic requirement is to 
make the estimation error “as small as possible”. 


desired 

response 


d(n) 


estimation 





error 

input 


output + 

km 

y(n) 

Linear discrete-time filter 

V(n) 


w - [ W 0 w ( W 2 ] 



Figure 2.1. Structure for the Wiener filtering problem. 


Generally speaking, the filter design can be optimized in the sense of minimizing 
a cost function , or index of performance. Among others, the possible choices are the 
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mean-square value, the expected absolute value, or the third (or higher) power of 
the absolute value of the estimation error [7, 17]. The choice of the minimum mean- 
square-error (MSE) criterion has attracted more research attention than the others, 
because it has the advantage of leading to tractable mathematics [7]. In particular, 
the MSE criterion results in a second-order dependence for the cost function on the 
unknown coefficients w in the impulse response of the filter. Moreover, the cost func- 
tion has a distinct minimum that uniquely defines the statistically optimum design 
of the filter. 

Assume that the filter input and the desired response are single realizations of 
jointly wide-sense stationary stochastic processes, and denote the FIR coefficient vec- 
tor of the filter as w = [w 0 wi ■■■ w M - i] T - The filter output u(n) and the estimation 
error e(n) at discrete time n are defined as 

M-l 

v(n) = w m y(n - m) = w T y(n), 

m=0 

and e(n) = d(n) - v(n) (2.13) 

where y(n) = [j/(n) y(n - 1) • • • y(n - M + 1)] T . To optimize the filter design, we 
choose to minimize the mean-square value of the estimation error, J = E[e(n)e*{n)\ = 
£[|e(n)| 2 ]. For the cost function J(w) to attain its minimum value, all the elements 
of the complex gradient vector VJ must be simultaneously equal to zero. (For a 
detail treatment of the concept of a complex gradient operator, please refer to Kay 
[23], Haykin [7], and Brandwood [24].) Under this condition, the filter is said to be 
optimum in the mean-square- error sense. 

Let R denote the M-by-M correlation matrix of the tap input vector y(n) = 
[ y (n) y(n — 1) * • • y{n — M + 1)] T in the transversal (FIR) filter of Figure 2.1: R = 
E[y(n)y H (n)}. Correspondingly, let P denote the M- by-1 cross-correlation vector 
between the tap input vector and the desired response d(n): P = E[y(n)d*(n)]. Let 
v 0 (n) denote the output produced by the filter optimized in the MSE sense, with e 0 (n) 
as the corresponding estimation error. The essential idea of Wiener filtering hinges 
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on two important results. First, the necessary and sufficient condition for minimizing 
the cost function J is that the corresponding value of the estimation error e 0 (n) is 
orthogonal to each input sample y(n) contributing to the estimation of the desired 
response at time n. Furthermore, it can be shown that under the same optimum 
condition, e 0 (n) and the estimate of the desired response v 0 (n) are orthogonal to each 
other. These statements constitute the principle of orthogonality and its corollary [7], 
which in mathematical terms are given as 

E[y(n - m)e*(n)] = 0, m = 0, 1, 2, • • • , M — 1, 
and £[u 0 (n)e’(n)] = 0. (2-14) 

Secondly, the optimum tap-weight vector (impulse response) of the transversal filter, 
denoted as w 0 = [u;^ w 0 i ••• w oM -i] t , can be obtained by solving the Wiener-Hopf 
equations 

Rw; = P. (2.15) 

Therefore, w* = R -1 P, assuming the correlation matrix R is nonsingular. 

2.4.2 LMS Algorithm 

When the filter’s performance index J = i?[|e(n)| 2 ] is a known function of w, 
Newton’s search method can be applied to minimize the required number of iterations. 
However, in many practical adaptive system applications the cost function J(-) is 
unknown and must be measured or estimated on the basis of stochastic input data. 
Among others, the method of steepest descent, which adjusts the filter weight vector 
in the direction of the gradient VJ at each iteration step, has thus far proven to be 
the most widely applicable. If it were possible to make exact measurements of the 
gradient vector VJ(n) at each iteration, and if the step size p is suitably chosen, then 
the tap-weight vector computed by using the steepest-descent algorithm would indeed 
converge to the optimum Wiener solution. In realty, however, exact measurements 
of the gradient vector are not possible, since this would require prior knowledge of 
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both the correlation matrix R of the tap inputs and the cross-correlation vector P 
[7]. Consequently, the gradient vector of J = £’[|e(ra)| 2 ] must be estimated from the 
available data. To simplify this problem, Widrow et al. [17, 25] use |e(n)| 2 itself 
as an estimate of J, and suggest that the required direction in which the weight 
vector should be changed is the opposite (or negative) direction of the maximum 
instantaneous rate of increase of the error power |e(n)| 2 with the weight vector. They 
come up with the LMS algorithm, which updates the filter coefficients according to 

e(n) = d(n) - w r (n)y(n) 

w(n-|-l) = w(n) -|- /xe(n)y*(n) (2.16) 

where w(n) = [u; 0 (n) uq(n) ••• u>M-i(rc)] r is the filter coefficients estimate at time 
instant n and p the adaptation gain or step size. 

In the LMS algorithm, the correction term /ze(n)y*(n) applied to the tap-weight 
vector w(n) at time n + 1 is directly proportional to the tap-input vector y(n). There- 
fore, when y(n) is large, the LMS algorithm experiences a gradient noise amplification 
problem [7]. To overcome this problem, we may use the normalized least-mean-square 
(NLMS) algorithm [26], a specific form of the LMS algorithm with a reparameterized 
step size, viz., 

e(n) = d(n) - w T (n)y(n) 

w(n + l) = w(n) + || y ^| --e(n)y , (n) (2.17) 

where p is the real positive step size. Using the idea of the projection algorithm [27] 
in the control literature, Slock [28] indicates that the NLMS algorithm is a potentially 
faster converging algorithm compared to the LMS algorithm, when the design of the 
adaptive filter is based on the usually quite limited knowledge of its input statistics. 
This performance advantage of the NLMS algorithm over its LMS precedent is also 
observed in the area of adaptive radar signal processing [29]. 
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Although the applications of adaptive filtering are quite different in nature, nev- 
ertheless, they have one basic common feature, i.e., an input vector and a desired 
response are used to compute an estimation error, which is in turn employed to con- 
trol the values of a set of adjustable filter coefficients. Dependent on the manner 
in which the desired response is extracted, the functions of the four basic classes of 
adaptive filtering applications can be categorized as: (1) system identification; (2) 
inverse modeling; (3) linear prediction; (4) interference (noise) canceling [7, 17]. 

Figure 2.2 shows the block diagram of a system identification configuration where 
an adaptive filter is used to provide a linear model that represents the best fit (in 
some sense) to an unknown plant characterized by its impulse response h. The plant 
and the adaptive filter are driven by the same input. The plant output supplies the 
desired response for the adaptive filter. If the plant is dynamic in nature, the model 
will be time- varying. Since the 1970s, many published papers have contributed to 
the understanding and confirmation of the LMS algorithm’s performance in tracking 
an unknown system, especially in the nonstationary signal environment [30, 31, 32]. 
In Chapter 4, through extensive simulations, we will evaluate the capability of adap- 
tive noise canceler (ANC) for the rejection of sinusoidal interference, following the 
configuration in Figure 2.2. Performance comparison between the LMS and NLMS 
algorithms will also be presented. 
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Figure 2.2. System identification via LMS adaptive algorithm. 
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CHAPTER 3 


HIGH-RESOLUTION FREQUENCY ESTIMATION VIA EM 

The essential ideas underlying the expectation maximization (EM) algorithm 
have been presented in special cases by many authors [33, 34, 35]. Dempster, Laird 
and Rubin [36] first recognized the expectation step (E-step) and the maximization 
step (M-step) in their general forms, and introduced it for computing maximum 
likelihood estimates from incomplete data. Since the late 1980s, the EM algorithm has 
attracted signal processors’ attention, and its significant contributions particularly to 
the area of sensor-array signal processing have been found in the literature [37, 38, 39]. 

In light of the frequency estimation problems addressed in the previous chapter, 
we develop a computationally efficient scheme for joint ML estimation of the signal 
spectral parameters, based on the iterative EM algorithm. This chapter begins with 
an overview of the data model and the basic ideas behind the EM algorithm, and then 
shows how this algorithm can be utilized to implement the ML frequency estimator. 
The important problem of determining the number of signals (order selection) is also 
addressed. To achieve the goal of correct order selection, we come up with an order 
recursive combined signal detection and estimation scheme via the EM algorithm. Fi- 
nally, through intensive simulation, we show how this algorithm attains high spectral 
resolution capability, low SNR threshold 1 , and high correct detection probability, a 
set of performance indices required of most modern signal processors. 


*SNR threshold is the lowest SNR level above which a frequency estimator approximately attains 
the performance of an ML estimator. 



3.1 Signal Model for the EM Algorithm 


3.1.1 Complete and Incomplete Data 

The EM algorithm, developed in [36], is a general approach to iterative com- 
putation of maximum-likelihood estimates when the observed signal samples can be 
viewed as incomplete data. The term “incomplete data” in its general form implies 
the existence of two sample spaces y and X and a many-to-one mapping x — ► y(x) 
from X to y. Instead of observing the “complete data” x in X, we observe the “in- 
complete data” y in y. Let the density function of x be / X (x; 9) with parameters 
9 6 f! and let the density function of y be given by 

fy{y\°)= j /x(x; 9)dx (3.1) 

JX(y) 

where X(y) — {x: y(x) = y}. 

3.1.2 Definition of EM Algorithm 

Given the observations Y = y, the MLE of 9 can be obtained by maximizing 
f y (y‘,9). However, in many statistical problems, maximization of the complete-data 
specification /r(x; 9) is simpler than that of the incomplete-data specification / y (y; 9). 
Following Dempster et al. [36], take logarithms of each side of 

/»( y; 9) = /*(x; 9)/f x |„(x; %) 

and then take conditional expectations given Y = y, under a parameter value 9k, to 
obtain 

L{y\9) = ln/ v (y; 9) 

= E x [lnf x (x;9)\Y = y,9 k ] - E xly [\nf x{Y=y (x;9)\Y = y,9 k ] 

= Q(9\9 k ) - H(9\9 k ). (3.2) 

It is a well-known consequence of Jensen’s inequality that H(9\9 k ) < H(9 k \9k) [36, 40], 
and thus implies 

ye Q(9\9 k ) > Q{9 k \9k), L(y\ 9) > L{y, 9 k ). (3.3) 


18 



Letting 0* denote the kth guess of the MLE of 0, we then have this iterative EM 
algorithm. 

Expectation (E-step): Determine the average log-likelihood of the complete 

data 

Q(0 \9 k ) = £ r [ln/*(x;0)|Y = y;0*] 

= J ln/ x (x; 0)/ x | y (x|y; 0jt)dx. (3.4) 

Maximization (M-step): Maximize the average log-likelihood of the complete 

data 

0 fc+ i = argmax<5(0|0jt). (3.5) 

$ 

At convergence we hopefully will have the MLE. This issue will be discussed later in 
Section 3.4. 


3.2 ML Frequency Estimator via EM Algorithm 

3.2.1 A Model for Signal Decompositions 
Consider again the observed signal y(n) in Section 2.1 which consists of p complex 
sinusoids corrupted by complex white Gaussian noise 

y(n) = £ A, exp0'2x/,n) + u>(n), n = 0, 1, • • • , N - 1 (3.6) 

» = 1 

where A,- = |Aj|e J ^‘ is the complex amplitude of the ith sinusoidal component, and 
w(n) is complex white Gaussian noise with zero mean and variance cr 2 . In the EM 
formulation, we call {y(n)j incomplete data. Suppose there exists the unobservable 
complete data {^(n), X2( n )> ‘ ‘ * i x p( n )} 

x,(n) = Ai exp(;27r/,n) + u,(n), n = 0, 1, • • • , N - 1; t = 1 • • -p (3.7) 

where v,(n) is complex white Gaussian noise with zero mean and variance cr,- 2 , and 
u,(n) is independent of Vj(n) whenever i ^ j. Let {y(n)} and (xi(n), x 2 (n), • • • , x p (n)} 
be related through the following noninvertible many-to-one mapping 

y( n ) = £ x »( n ) » ( 3 - 8 ) 

«=i 
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which implies 


p p 

w(n) = ^2 v i(n) and <r 2 = y cr, 2 , 

i=i <=i 

P P a ? 

or equivalently, y ft, = y — jp = 1. (3.9) 

i=i »=i a 

3.2.2 Frequency Estimation via EM: a Brief Derivation 
Let x t - = [x,(0) x,(l) • ■ ■ Xi (N- 1)] T , e, = e, •(/,•) = [1 e> 2 **' 2 • • • e ;Wd*-i)]T 

Oi = { fi,Ai }, X = [xi x 2 ••• x p ], 0 = {0 u 6 2 ,---6 v }, and 0* = {Q lk , 0 2 „, • • • 9 Pk ), the 
itth guess of the MLE of 0. From (3.7) we have, upon noting the independent data 
set assumption, 

In /r(X; 0) = tlnfiXiA) 

1 = 1 
P 

= 5 > 

1 = 1 

= l*i( n ) “ ^.exp(j27r/,n)| 2 

1=1 n =0 

= 9-^2 ~2 llx,- - A,-ei|| 2 (3.10) 

i=l C i 

where g is a constant independent of the parameter set 0. 

Given the observed signal y and previous parameter estimates 0*, by taking the 
expectation of (3.10) we find 

Q(0|0*) = £*[ln/ r (X;0)|y;0 fc ] 

= E(g\y; Q k ) - £ l|£(x,-|y; Ok) - ^e t -|| 2 

1=1 a i 

= E{g\ y;e*)-£-WWi) (3.11) 

.=1 a i 

where 0{ k = {fi k ,Ai k } is the Hh MLE guess of 0,. Using the standard result for 
conditional expectations of jointly Gaussian random vectors, and following the similar 
derivation by Kay [23] for the case of real sinusoidal signal, it can be shown that 

x, = E(x,|y;0)t) 

= (3.12) 
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where e tJk = e, t (/,J = [1 e j2,r/, <=e- ?2,r/ ’* 2 ••• e - ,2,r -M Ar_1) ] T , and x,- = [x,(0)x,(l) ••• x,(fV — 
1)] T . Note that £(x,|y; 0*) can be thought of as an estimate of x,(n) in the complete 
data set since from (3.12), 

Xi(n) = A ik exp (j2irf ik n) + ^ ^y(n) - A ik exp(j2irf ik n)^j . (3.13) 

Obviously, in (3.11), to maximize <5(0|Ojfc) with respect to 0 is to minimize each 
S{(Ai,fi) individually, knowing that E{g\y\ 0t) is independent of 0. According to 
Section 2.2.1, minimization of the scoring function for single sinusoidal 

parameter estimation, can be achieved by choosing 

. 1 [ IT ^ |2 

/.•*+! = ar g max — | e , x,| 

and A ik+i = -^e£ +i x,-, for i = 1, • • • ,p. (3.14) 

Up to this point, it can be clearly recognized that (3.13) (E-step) and (3.14) (M-step) 
constitute the iterative EM algorithm for ML frequency estimation. 

Notice that (3.14) represents a one-dimensional MLE processor producing its by- 
product Ai k+i as the ML amplitude estimate of the single- sinusoid in x,(n). According 
to Section 2.1.2, given the current frequency estimate f* +1 = [fi k+1 h k+l ■ • • f Pk+1 ], 

A t+1 = (E" .E^.r'E^.y 

is the amplitude estimate that will maximize the joint log-likelihood L(y; fit+i), where 
Efc+i = [e ljfc+1 e 2 fc+1 • • * e Pfc+ J. In other words, using this amplitude estimate instead 
of that produced by (3.14) will lead to a more generalized EM (GEM) algorithm as 
described in Dempster et al. [36]. 

Finally, since the erf ' s are not unique, they can be chosen arbitrarily as long as 
(3.9) is satisfied. Recognizing that after the &th iteration, i ]; k = is the SNR 

estimate of the expected ith sinusoidal component x,(n), it is quite reasonable to 
choose erf such that 

Tft k = rj 2 k = • • • = rj Pk = constant 
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for each MLE processor in (3.14). Considering this choice of of in (3.9), it is very 
straightforward to show 


Pi = rk = 


\Ai 


r, for * = 1,2, 


n=iiAj 

Now, the generalized EM frequency estimation algorithm is summarized below. 

• Initialization: Given the initial frequency estimate f 0 = [/i 0 f 2o • • • / Po ], and thus 
the initial amplitude estimate A 0 = (E^E 0 ) _ 1 E^y, continue the following EM iter- 
ation until ||f*+i - ffcHoo < e, where ||h||oo = max,- 

• Expectation (E-step): Decompose y(n) into a set of expected sinusoidal compo- 
nents (x,(n), i = 1 • • • p} 

l*J 


x,(n) = A ik exp(j2ir fi k n) + 


or x 


TLi IA J 

' = + eL’ki ( y " S' 4< ‘ e< ‘) ' 


y( n ) - Ai * exp(j27r/, fc n) , 


1 = 1 


(3.15) 


• Maximization (M-step): Maximize the log-likelihood of each expected sinusoidal 
component separately by finding 


1 I i2 

/.•*+, = ar g max -y , * = 1, 2, • • • , p 

and Ajt+i = (Ef +1 E fc+ i) -1 Ef +1 y. 


(3.16) 

(3.17) 


At its convergence, denote the final frequency and amplitude estimates as f p = 
[fi h ' ' - fp] an( i A p = [Ai A 2 • • ■ A P } T . Define the pth-order model residual x p as 

z p = y -^A,e. (3.18) 

1=1 

where e, = e,(/,) = [1 e* 2 *** e j2lt ^ 2 e- ,2ir ^d Ar - 1 )] T . From now on, for the sake 
of notational clarity, we will refer to this EM frequency estimator by the following 
formula: 

{f p , Ap,z p } = £M(y,p, {f 0 ,A 0 }). 
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3.3 Combined Signal Detection and Estimation 

Like most iterative frequency estimation algorithms, the EM algorithm requires 
a reasonably accurate initial parameter estimate 0o to begin with. Taking a closer 
look at (3.15), (3.16) and (3.18), it is not difficult to observe some inherent “order- 
recursive” features associated with this algorithm. Let p be the true model order 
(number of sinusoidal signals), and m be the order we pick. First notice that when 
m = 1, Xi = y according to (3.15), and /i is the frequency at which the periodogram 
of y attains its maximum. Secondly, the model residual z m defined by (3.18) reveals 
to us the signal’s structural remnant yet to be modeled, when m < p. Based upon the 
above observations, we present the order-recursive EM algorithm, for ML frequency 
estimation as follows: 

• Let m = 1, and x m = y. Find /j, A\ and Zi using (3.16)~(3.18). fi = [fi\. 

• For m = 2 to p 
do 

1. initialization: 

| I JJ I 2 a 

/m 0 = arg max — e"z m _J . Let f 0 = [f m _i f mo ]. 

Jm iV ' 1 

Calculate Ao = (Ej^Eo) _1 E^y, where Eo = Eo(fo). 

2. EM iteration: 

{ f m i A m > 2 m ) m, {fi), Ao })• 

end. 

At this point, one can see that an interesting feature this algorithm possesses is its 
built-in initialization procedure. For each model order m < p, the initial frequency 
estimate is the optimal (in the ML sense) estimate of previous order (m — 1) plus one 
additional estimate extracted from the model residual z m -i- 

When the true model order p is unknown, some sort of “stopping rule” has to be 
applied to the order-recursive EM algorithm. Apart from the information theoretic 
criteria mentioned in Section 2.3 for order selection, a natural way to stop modeling 
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the sinusoids-plus-noise time series is to test the whiteness of the model residual 

z m . Since order increment and frequency estimation in this algorithm hinge on the 

periodogram of the model residual, periodogram-based methods might be used to 

construct statistical tests for departure from white noise. If the suspected source of 

departure is the presence of a single sinusoidal component at unknown frequency /, a 

natural test statistic is the maximum periodogram ordinate [41], T = max/ ^-|e w y| 2 . 

Let {/fc} be the periodogram ordinates, where Ik = ^|e^y| 2 for k = 1,2, -,7V. 

Under the hypothesis that y is white noise, the distribution function of T requires 

the knowledge of the noise variance a 2 . However, in most real world applications, cr 2 

is unknown. To solve this problem, Fisher [42] devised the exact distribution function 

of To = , _ jv , the ratio of the maximum to average periodogram ordinate. In 

77Z_,*=i 7 * 

particular, he showed that 

R f AM 1 

P(T„ > Ng) = g [ t , (A ,_ i;)! ] (-1)‘-'(1 - t?)"- 1 (3.19) 

where R is the largest integer less than q~ l . Based upon the above arguments, a 
reasonable way to decide whether to increase the model order is to test the whiteness 
of model residual z m using (3.19). In other words, if the model order we pick is correct, 
the residual series z p is approximately a partial realization of white noise, and thus 
passes the whiteness test. Figure 3.1(a) shows the block diagram of the EM frequency 
estimator for a given order p, and Figure 3.1(b) provides a flow-chart representation 
of the order-recursive EM frequency estimator with Fisher’s To statistic as the model 
order selector. The comparison of order decision performances for Fisher’s To statistic 
and information theoretic criteria like AIC and MDL will be given in the next section. 

3.4 EM Convergence Properties and Simulation Results 
In general, if the log-likelihood function L(Q) (i.e., L( y; 0)) has several (local 
or global) maxima and stationary points, convergence of the EM sequence {0jt} to 
either type of point depends on the choice of starting points. This phenomenon has 
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(a) 



(b) 


Figure 3.1. (a) Block diagram of EM frequency estimator, and (b) flow chart repre- 
sentation of the combined order-recursive frequency estimation and signal detection 
algorithm. 
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been well recognized by statisticians and signal processors [36, 37, 38, 43]. The most 
striking characteristic of the EM algorithm is that L(Qk+ i) > L(Ojt) for any EM 
sequence {0*}. This implies that if L(Qk) is bounded above, then it converges to 
some L*. Miller and Fuhrmann [38] suggest that in real world applications, L(Qk) 
can be upper bounded by properly choosing a closed and bounded parameter space Cl 
for 0. More specifically, Wu [43] shows that, if (3(0|$) in (3.11) is continuous in both 
0 and $, a condition satisfied in most practical situations, then all the limit points 
of any instance {0*} of an EM algorithm are stationary points of L. Furthermore, 
1,(0*) converges monotonically to L * = L(0*) for some stationary point 0*. For a 
rigorous treatment of EM convergence properties, please refer to Wu [43] and Boyles 
[44]. 

In this section, we demonstrate via simulation that the EM algorithm for ML 
frequency estimation will resolve signal components in situations of small sample 
size and low SNR which cause other high resolution estimators to fail. We consider 
two signal scenarios that consist of two or three sinusoids with different frequency 
separation and SNR level. For each SNR level of interest, 200 Monte Carlo simulations 
are undertaken. In the first case, the observed signal is generated as 

y(n) = \A\e** hn + \A\e jQ - 2Sx e j2 * hn + w(n), n = 0, 1, • • • , 24 

where |A| is the scalar amplitude, /i = 0.19, /j = 0.21, and w(n) is a complex white 
Gaussian sequence. To help visualize the convergence behavior of EM iterations, 
Figures 3.2 and 3.3 show the log-likelihood surfaces, surface contours, and trajectories 
of EM frequency estimation starting from different initial estimates. It is the result 
of processing a snapshot of y(n ) at both low SNR (0 dB) and high SNR (20 dB) 
individually. 

Similarly, in the second case, the observed signal is 
y(n) = \A\e’ 2 * hn + |A|e jtU5 V 2ir/sn + |A|e j0 - 85 V 2,r/3n + tn(n), n = 0, 1, • • • , 24 
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log-likelihood surface & trajectories of EM iterations 



Figure 3.2. Log-likelihood surface, surface contours, and trajectories of EM iterations 
with different initial frequency estimates and f\ = 0.19, / 2 = 0.21, SNR = 0 dB, 
N = 25, EM estimate : /j = 0.1874, / 2 = 0.2136. 






log-likelihood surface & trajectories of EM iterations 



Figure 3.3. Log-likelihood surface, surface contours, and trajectories of EM iterations 
with different initial frequency estimates and fi = 0.19, f 2 = 0.21, SNR = 20 dB, 
N = 25, EM estimate : fi = 0.1909, f 2 = 0.2100. 





where f 3 = 0.365, a frequency further away from the first two. Notice that with 
a sample length of 25, the Fourier resolution bandwidth is 1/25 = 0.04, which is 
larger than the frequency spacing of the first two sinusoids. Also, none of the se- 
lected signal frequencies is a multiple of this Fourier resolution bandwidth. In terms 
of being unbiased, small mean square error (MSE), and low SNR threshold (about 
2 dB), Figures 3.4 and 3.5 reveal the excellent performance of our algorithm, where 
the Cramer- Rao bound [3, 45] is used as the benchmark for performance assessment. 
Although global convergence is still an issue, we show through simulation that when 
the initial frequency estimate of each component is within approximately one resolu- 
tion bandwidth of of the global maximum, convergence can be achieved. In fact, the 
optimal built-in initialization procedure of our algorithm, as pointed out in Section 
3.3, has successfully carried out this task. 

Finally, we demonstrate via simulation the modeling (detection) performances of 
Fisher’s T 0 statistic, MDL and AIC, operating in conjunction with the EM frequency 
estimation algprithm. The simulated signal is the f/tree-sinusoids-plus-noise described 
above. For the information theoretic criteria (MDL and AIC), model order selected 
for examination ranges from 1 to 5. Therefore, we count orders 1 and 2 picked by 
these criteria as underfitting, and orders 4 and 5 as overfitting. In the use of Fisher’s 
T 0 statistic, H 0 (noise hypothesis) is rejected at the significance level a = 0.01, which 
equivalently thresholds T 0 at 6.9547 according to (3.19) with N = 25. Figure 3.6 
shows the model selection capabilities of these criteria, based upon 100 Monte Carlo 
simulations at each SNR level of 0 dB, 3 dB, 5 dB and 10 dB. Notice that at all SNR 
levels of interest, MDL consistently outperforms Fisher’s T 0 statistic with higher 
probability of correct modeling, and AIC completely overfits the signal model in all 
simulations. 
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Figure 3.5. Performance of the EM algorithm for ML frequency estimation, p 
fi = 0.19, f 2 = 0.21, f 3 = 0.365, N = 25. 
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Figure 3.6. Comparison of detection performances for Fisher’s To statistic (level of 
significance a = 0.01), MDL, and AIC. 
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CHAPTER 4 


LMS ADAPTIVE FILTER FOR SINUSOIDAL PROCESS 

The suppression of a sinusoidal interference corrupting an information-bearing 
signal is a problem often encountered in many signal processing applications. The 
traditional way of dealing with this problem is to design a fixed notch filter tuned 
to the frequency of the interference. To design the filter, precise knowledge of the 
interfering signal’s frequency is always required. When the notch is desired to be 
very sharp and the sinusoidal interference is known to drift slowly, a fixed filtering 
approach may have difficulty solving the problem. It has been shown by Widrow et 
al. [46] and Glover [47] that a notch filter realized by an adaptive noise canceler can 
offer advantages such as easy control of notch bandwidth, an infinite null, and the 
capability of adaptively tracking the exact frequency and phase of the interference. 

In this chapter, we will investigate the performance of the adaptive noise canceler 
in the suppression of multiple complex sinusoidal interferences. Different functions 
of the adaptive noise canceler, e.g., notch filter, a process decorrelator and a line 
enhancer, as well as the convergence properties associated with various model orders 
of both the signal and system will also be addressed. 

4.1 Notch Filter Realized by Adaptive Noise Canceler 

Figure 4.1(a) shows the block diagram of a dual-input adaptive noise canceler 
(ANC). The primary input supplies an information-bearing signal and an interfering 
noise of multiple sinusoids that are uncorrelated from each other. The reference input 
consists of a correlated version of the sinusoidal interferences. For the adaptive filter, 
we use a transversal filter whose weights are adapted by means of the LMS algorithm. 
The reference input and the primary input are given respectively as 

p 

V( n ) = exp(ju>,n), 

1 = 1 


(4.1) 



Primary input System output 



(a) 



Figure 4.1. (a) Schematic representation of adaptive noise canceler and (b) equivalent 
model in the z-domain. 
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(4.2) 


P 

and d(n) = s(n) + '^ / Bi exp(ju t n) n = 0, 1,2, ••• 

i=i 

where s(n) is the information-bearing signal, A , and B{ are complex amplitudes, and 

u>i = 2irfi (0<fi< 1). 

The filter will use the reference input to provide (at its output) an estimate of the 
sinusoidal interfering signal contained in the primary input. Thus, by subtracting the 
adaptive filter output, v(n), from the primary input d(n), the effect of the sinusoidal 
interference is diminished. According to Section 2.4.2, the LMS algorithm updates 
the tap-weight as follows: 

M-l 

v(n ) = w m(n)y{n - m) (4.3) 

m — 0 

e(n) = d(n) — v(n) (4-4) 

w m (n + 1) = w m (n) + ny‘(n - m)e(n), m = 0, 1, • • • , M - 1 (4.5) 

where M is the total number of tap weights (order of filter) in the transversal filter, 
and \i the constant step-size parameter. Let V(z) and E(z ) denote the z-transform 
of the filter output u(n) and the estimation error e(n), respectively. Following the 
schematic representation in Figure 4.1(a), we may lump the sinusoidal reference input 
y(n), the transversal filter, and the weight-update equation of the LMS algorithm into 
an open-loop system defined by a transfer function G(z) = ^fj", as in the equivalent 
model of Figure 4.1(b). Our goal is to find F(z), and thus G{z), given E(z). Starting 
from the weight-updating formula in (4.5), by taking the z-transform of both sides, 
we get 

zW m (z) = W m (z) + nZ{y m {n- m)e(n)}, m = 0, 1, • • • , M - 1 (4.6) 

where W m (z ) is the z-transform of u; m (n). Substituting (4.1) into (4.6), we find 

W m (z) = [£ Ale^E(z^)] (4.7) 

1 z U=i 

where E(z) is the z-transform of e(n), and thus E(z e JWi ) is E(z) rotated clockwise 
around the unit circle through the angle u>j. Furthermore, taking the z-transform of 


35 



u(n) in (4.3), with y(n) and W m (z) replaced by (4.1) and (4.3) respectively, it can be 
shown that 


V{z) = Z\ E u )m (n)^A l e^' m e^ n 


M - 1 p 


= 5;Dv"7{».(»r) 


m=0 i=l 


M— 1 P 
m=0 t= 1 


M - 1 p 


* » //?r ^ P / ■/ \ 

= E DAe'*") i 


_ V- KM I 

k 1 - 2 - 


/lAflAipz-V" 4 


+ EE 


TI component 

^ y- l c* u,% g, ( ze i{u k -w t )'' 

s 2 -j i ^ z -i e ju >i \ > 


i=l k = \ 
k*i 


TV component 


where 


ivt — i. in-1 

Pik{M) = E e j(w *- u ' ,)m = El e j2 * Uk ~ f ' )m 

m=0 m=0 

= Sm(( A~ /, ) M r eJfA - /,)(Af - I) ^ (4.9) 

sm((/jt - /,» 

Taking a closer look at the expression for V(z) in (4.8), we consider the first 
term as a time-invariant (TI) component , and the second term as a time-varying 
(TV) component [47]. According to (4.8), the effect of the time- varying component 
depends on the factor 0ik(M) defined by (4.9). Particularly, when « 0, G(z) 

is determined by retaining only the time-invariant component of V(z ), and the ANC 
behaves like a fixed multiple- notch filter. Letting 6/ = min/t^, 1 fk ~ fi\i the minimum 
frequency spacing between the interfering sinusoids in the reference input, the time- 
invariance condition of the ANC can be satisfied by choosing M > j-. The open-loop 


transfer function G(z) is therefore 


G(z) = 


E(z) ti 


k i 


p.M\Ai\ 2 z 1 e : 


z~ l e }Wi 


(4.10) 
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and the notch filter realized by Figure 4.1(a) can be recognized as a closed-loop feed- 
back system with transfer function 


H(z) = 


E(z) 1 

D(z) 1 + G(z) 


1 

1 + ^2 ^\ A >\ 2z ~ leM 

i=i 


(4.11) 


1 — z~ x e }w ' 

From Equation (4.11), the zeros (notches) of H(z ) are at the poles of G(z); that 
is, they are located on the unit circle at e■ ?u ' , . Furthermore, if a small value of the 
step-size parameter p is chosen, such that pM\Ai\ 2 < 1, i = 1,2, ■ • • ,p, the poles of 
H(z) can be approximately located at 


Zi « (1 — pM\Ai\ 2 )e- 




This fact implies that the poles of H(z) lie inside the unit circle, and thus that the 
ANC is stable, as it should be for practical use in real time [7]. 

Following Equation (4.11), Figures 4.2(a) and 4.2(b) give the time-invariant por- 
tion of an adaptive notch filter’s frequency response and its corresponding pole-zero 
plot, with | Ax | = \A 2 \ = \A 3 \ = 1, [/, / 2 / 3 ] = ( 0.38 0.45 0 . 76 ], M = 32 , and a 
very small adaptation rate p = 0.0002 to narrow the notch bandwidths. Due to this 
choice of pM\A{\ 2 , as can be seen in Figure 4.2(b), poles and zeros of H{z) tend to 
overlap each other on the unit circle. Furthermore, Figure 4.2(c) provides the result 
of an experiment performed to characterize the adaptive notch filter’s response to a 
complex unit impulse function as the primary input, i.e., d{n) = ^(1 + j)S(n). In 
Figures 4.2(a) and 4.2(c), the frequency response of H(z) and the spectrum of e(n ) 
are evaluated at 500 normalized digital frequencies from 0 to 1. 

4.2 Adaptive Noise Canceler as Process Decorrelator 
In Section 4.1, we assume that the reference signal is deterministic while trying 
to implement a multiple-notch filter via ANC. However, to deal with signals encoun- 
tered in real world applications, it is more often the random nature associated with 
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Tl frequency response of adaptive notch filter 



pole-zero plot of H(z) 



Spectral magnitude of e(n), given d(n) - unit impulse function 



Figure 4.2. (a) Time- invariant frequency response of a notch filter realized by ANC; 
(b) pole-zero plot of H(z)- (c) spectrum of e(n) when d(n) is a unit impulse function. 
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the signals that necessitates adaptive approaches. Assuming that all the signals con- 
sidered are realizations of some random processes, and only limited knowledge of the 
correlation properties among the signals is available, we can use the ANC as a process 
decorrelator to extract the information-bearing signal. Let the reference signal y(n ) 
in Figure 4.1(a) be represented as 

p 

y(n) = '%2A i exp(j2Trf i n + <t>i) + a(n), n = 0, 1, • • • , N - 1 (4.12) 

1=1 

where A{ are constant amplitudes, phases <f>{ are independent random variables uni- 
formly distributed over [0 2 tt), and a(n) is an independent CWGN with zero mean and 
variance It has been shown that y(n) is a wide-sense stationary (WSS) process 
with autocorrelation function (ACF) given as [3] 

r yy (k) = £ A 1 exp(j2Trfik) + <T 2 a 6(k). (4.13) 

1=1 

Denote S as the diagonal matrix with the power of the ith sinusoid, Si = Af, as its 
ith diagonal element so that the M x M autocorrelation matrix for y(n) is 

Ryy = £ Si**" + °l l = ESE " + (4-14) 

1=1 

where E = [ei e 2 • • • e p ] with e, = [1 e j2 ’ r ^ e j2ir *' 2 ■ • • e j2 ’ r -f^ M-1 )] T . Notice also that 
Ryj, is Hermitian. 

Suppose that the primary input d(n ) and the reference input y(n ) are related as 
follows: 


d(n ) = s(n) + x(n) 

L- 1 

x ( n ) = ^2 h k y(n- k) = h T y L (n) (4.15) 

k = 0 

where h = [/io hi • • • hi-\] T characterizes the linear correlation between x(n) and 
y(n), s(n ) is an information- bearing signal process uncorrelated with x(n), and yz,(n) = 
[j/(n) y(n — 1) • • • y(n — L + l)] r . We call L the order of correlation between the 
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primary input and the reference signal. Letting P = £[yM(ra)d*(n)], the M x 1 cross- 
correlation vector between the tap input vector and the primary input, and based 
upon the above assumptions, it is straightforward to show that 

P = Rjv/Lh* 

where R ml = E[y M{n)yi (n)]. Therefore, assuming that the autocorrelation matrix 
Rj, y is nonsingular, the optimum Wiener solution described by (2.15) can be obtained 
as 

W 0 = RyyRM£,h*. (4-16) 

Up to this point, it is quite obvious that when M = L, Ryj, = Ra/l, and thus w Q = h. 

To show the performance of the ANC functioning as a process decorrelator, 
we undertake another experiment of 1000 Monte Carlo simulations. The reference 
signal is generated according to (4.12) with [/i / 3 ] = [0.38 0.56 0.75] and 10 dB of 

sinusoids-to-noise power ratio. For each realization of y(n ), d(n ) is correspondingly 
generated by (4.15), where s(n) is replaced by a white Gaussian noise z(n) with power 
30 dB below that of x(n). We choose the correlation vector h as 

_ j- g (-0.2+i0.60487r)0 e (-0.2+j0.6048»)l g (-0.2+j0.6048w)2jT 

= [1.0000 +j0 — 0.2647 + >0.7748 - 0.5302 -j0.4102] r , 

and assume that the number of filter taps is equal to the order of correlation, i.e., 
M = L = 3. If z(n ) is taken as the “plant noise”, we are equivalently facing a system 
identification problem as depicted in Figure 2.2. 

Figure 4.3 demonstrates the ensemble averages of both the filter’s tap weight 
vectors and the power of the estimation error e(n). At the end of all simulations, we 
obtain the filter’s tap weight vector estimate 

w 0 = [0.9979 - >0.0005 - 0.2679 + >0.7746 - 0.5324 - ;'0.4096] T , 
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ReaJ part of ANC FIR parameter vector, p - 3 



Imaginary part of ANC FIR parameter vector, L ■ 3 
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Figure 4.3. Process decorrelator via ANC, M = L — 3, /i = 0.0011: evolution of tap 
weight vector: Re{w(n)} and Im{w(n)}; learning curve: 10 log(J(n)) vs. time. 
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which is almost equal to w„ = h, the Wiener solution of this case. Under such an opti- 
mum condition, only the plant noise is expected to be left at the process decorrelator’s 
output, i.e., e(n) z(n). In other words, by examining the “learning curve” at the 
bottom of Figure 4.3, one can find that the mean-squared error sequence, E[\e(n)\ 2 \, 
converges to about —30 dB, a phenomenon in agreement with the assumption made 
for the primary input signal. 

4.3 Adaptive Filter Design Considerations 
In Section 4.2, we assume that the reference signal y(n) and its linearly correlated 
signal x(n) are random-phased sinusoids mixed with white Gaussian noise, as modeled 
by (4.12) and (4.15). Due to the presence of the noise’s autocorrelation matrix, a 2 w I, 
at the right-hand side of (4.14), Ryj, is of full rank and nonsingular. These facts 
justify the existence of a nontrivial solution to the Wiener-Hopf equation, as expressed 
by (4.16): w* = R^'R^lE*. Moreover, an adaptive noise canceler is expected to 
suppress the unwanted interference x(n) by adaptively tracking its tap weight vector 
w(n) to this optimum w 0 . As pointed out before, when M = L, we have w 0 = h. 
However, in most applications, the order of correlation L is unknown to the signal 
processors. In this section, we will discuss the issue regarding the choices of two 
important design parameters: the filter’s number of taps, M, and the LMS (NLMS) 
algorithm’s adaptation rate /i. 

4.3.1 Choice of Filter’s Number of Taps M 
According to Orfanidis [48], the rule of choosing M with respect to L is that 
the adaptive filter must have at least as many delays as that part of d(n ) which is 
correlated with j/(n), in other words, M > L. Taking a closer look at Equation (4.16) 
can help us see why it is so. First, consider the case of overmodeling , i.e., M > L, 
then partition Ry V and its inverse matrix R“ v * as 

Riw = [R-a/l B] and R^ 1 = 
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where B is an M x (A/ — L) matrix, R,„ an L x M matrix, and R,; an ( M — L) x M 
matrix. Since R“ y ‘R vy = I, it is straightforward to show that R, u Rjv/z, = I and 
Rt/Rml = o, thus the Wiener solution can be found as 


w* = R,„! R\/i,h* = 


yy 


Ilxl 

0(M-L)xL 


h* = 


h* 

0 


(4.17) 


where 0 is the (M — L) x 1 zero vector. Therefore, according to (4.3) ,(4.4), (4.15), 
and (4.17), if w(n) = w 0 = [h 0 h x • • • 0 • • • 0] r , 


L - 1 

e(n) = d(n ) — wjy(n) = d(n) — ^ hky{n — k) = d(n) — x(n ) = s(n). 

k=0 


This implies the complete cancelation of the j/(n)-dependent part of d(n). 

When the adaptive filter is undermodeling the correlation between y(n) and d(n), 
i.e., M < L, following (4.15) we can express d(n) as 


d(n) = hfyi(n) + h^y 2 (n) + s(n) 


where hi = [/i 0 h x ■ ■ ■ h M -i] T , h 2 = [h M h M +i ••• yi(n) = [y(n) y(n - 

1) • • • y{n - M + l)] r , and y 2 (n) = [j/(n - M) y(n - M - 1) • • • y(n - L + 1)] T . 
Furthermore, it can be shown that 


w; = R ‘RMih* = RnfRn R !2 


hi 


— hj + R n Ri2h 2 , 


(4.18) 


and the optimum estimate of d(n) given yi(n) is 


d(n) = F[d(n)|y!(n)] = F[d(n)yf(n)]R n 1 y 1 (n) = wjy t (n) 


where Rj, y = Rn = £'[yi(n)y( / ] and R J2 = £[yi(n)y^]. More specifically, as w(n) 
converges to w 0 , 

d(n) = v(n) = wfy i(n) 

= hfyi(n) + h|’R( / 2 (Ri' 1 1 ) // y 1 (n) = hfy^n) + h^^^y^n), 


and thus the estimated information-bearing signal 

e(n) = d(n) - d(n) = h[[y 2 (n) - y 2/1 (n)] + «(n) 
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where y 2 /i (n) = £(y 2 (»)yf(n)]£[yi(n)y"(n)]~ 1 yi(n) = R 2 iR u 1 yi(”) is the opti- 
mum estimate of y 2 (n) based on yi(n). This analysis shows that the yi(n) part 
is completely removed from the primary input, and the y 2 (n) part is suppressed as 
much as possible. 

4.3.2 Choice of Adaptation Rate p and Signal Statistics 

In designing adaptive filters via the LMS algorithm, a problem of many concerns 
is the convergence behavior in tracking the optimum Wiener solution, where the 
adaptation rate /x plays a very critical role. Before delving into this issue, we need to 
know that the LMS algorithm is an example of a multivariable nonlinear stochastic 
feedback system , and such combined presence of nonlinearity and randomness makes 
its convergence (stability) analysis a difficult mathematical task [7, 17]. To alleviate 
the mathematical intractability in the convergence analysis of the LMS algorithm, a 
set of fundamental assumptions needs to be followed: 

1. The tap-input vectors y(l),y(2), • • • ,y(n) are statistically independent. 

2. At time n, the tap-input vector y(n) and the desired response d(n ) are 
statistically dependent, but independent of their previous counterparts. 

3. y(n) and d(n) are jointly Gaussian- distributed random variables for all n. 

The statistical analysis of the LMS algorithm based on the fundamental assumptions 
is called the independence theory. Please refer to Gardner [49] and Haykin [7] for 
a detail account. Our discussion regarding the choice of /x and its corresponding 
mean-squared error J(n) = f?[|e(n)| 2 ] will rely on the independence theory. 

Letting A,-, i = 1,2, • • • , M, denote the eigenvalues of the correlation matrix Ry y , 
the mean-squared error J(n) converges to a steady-state value J( oo) if, and only if, 
the adaptation rate /x satisfies 

0 < /x < (4.19) 

^max 

and <4 - 20) 
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where M is the number of filter taps, and A max is the largest eigenvalue of R^. Since 
R yy is usually unavailable to the ANC designers, when p is small compared to , 

^nux 

(4.20) can be simplified to a rule of thumb as 


0 < < 


(4.21) 


total input power 

where the total input power is an estimate of £’[y i/ (n)y(n)] = trfRyy]. When these 
conditions are satisfied, the LMS algorithm is expected to converge in the mean- 
square sense. Furthermore, define the minimum mean-squared, error as the MSE 
produced by the optimum Wiener filter, i.e., J ^ = E[\d(n) - wjy(n)| 2 ]. According 
to the independence theory, the mean-squared error produced by the LMS algorithm 
has the final value 

d min 


J(o o) = 


1_£ ’ 


(4.22) 


X 2-/iA, 

which is always in excess of the minimum value Jnun due to the variance of w(n) with 
respect to the optimum Wiener solution w 0 . A quantitative measure of this cost of 
adaptability is the misadjustment M, defined as 


M 


M a J{ oo) - J n 


Jrr 


t=l 

M 

1-EaA,/(2-/xA.) 
1 = 1 


(4.23) 


In Section 4.3.1 we see the effect of choosing M in general Wiener filtering. 
So far as the LMS algorithm is concerned, Equations (4.20) and (4.23) show that 
three principal factors affect the convergence behavior of this algorithm: the step- 
size parameter p, the number of filter taps M, and the eigenvalue distribution of the 
correlation matrix R yy . The condition number of R^, defined as 

A m j, v 


X(Hyy) — 


V . ’ 

/ 'min 


is a good indicator of the signal statistics. This ratio is commonly referred to as the 
eigenvalue spread (ES) or the eigenvalue disparity , a factor that controls the LMS 
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algorithm’s convergence speed [7, 48]. A large eigenvalue spread in the correlation 
matrix R yy corresponds to a highly self-correlated signal y(n). While a reference signal 
y(n) of such nature tends to slow down the LMS algorithm’s convergence process, the 
tracking dynamics of the NLMS algorithm appear to be significantly less sensitive to 
a variety of input signal distribution aspects than holds for its precedent [28, 50]. This 
point can be partially validated by the principle that the NLMS algorithm converges 
if and only if 

0 < ft < 2 (4.24) 

which is a condition on the step size parameter that is independent of the signal 
statistics. The fastest convergence occurs for 

A = 1, 


corresponding to the projection interpretation discussed in [27]. In next section, we 
will demonstrate the effect of these design parameters, and the performance compar- 
ison for the LMS and NLMS algorithms through simulation. 

4.3.3 Simulation Results 

Clearly, there are many practical problems for which the reference input process 
and the desired response do not satisfy the fundamental assumptions. An example will 
be the signal model assumed in Section 4.2 for the process decorrelator. Nevertheless, 
experience with the LMS algorithm has shown that the independence theory retains 
sufficient information about the structure of the adaptive process for the results of 
the theory to serve as reliable design guidelines, even for some problems having highly 
dependent data samples [7, 51]. 

To demonstrate the effects of signal statistics, correlation modeling, and the 
choice of /x, we perform three experiments of 1000 Monte Carlo simulations. In each 
experiment, the reference input and the primary input are generated according to 
(4.12) and (4.15) with L — 4, N = 400, and «s(n) being white Gaussian noise of power 
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level 30 dB below that of the sinusoidal interference. The sinusoid-to-noise ratio at 
the reference input is 5 dB. Performances of both the LMS and NLMS algorithms 
are examined in every simulation. Moreover, we choose the filter’s number of taps 
as M = 3, M = 4, and M = 5 for undermodeling, correct-modeling , and overmodel- 
ing the correlation between both inputs, respectively. To expedite the performance 
comparison, step size ji of the NLMS algorithm is fixed at 1 for all cases. 

In the first experiment, we have h = [h 0 hi h 2 /i 3 ] t with h k = e~ ak where 
a — 0.2 — j0.60487t, and three interfering sinusoids at frequencies [/x f 2 / 3 ] = 
[0.23 0.45 0.78]. Figure 4.4 shows the LMS algorithm’s convergence behavior in terms 
of its time evolution of filter’s weight vectors and output mean-squared errors. Notice 
that in the undermodeling case ( M < L), w(n) very quickly converges to the Wiener 
solution, but J(n ) converges to a steady level of -10 dB due to the reason addressed 
in Section 4.3.1 and the inherent misadjustment associated with the LMS algorithm. 
However, as M increases, the convergences of w (n) and J(n) are significantly slowed 
down. When M > L, the observed J(oo) is about —28.55 dB, very close to —28.66 
dB as predicted by (4.22). Similar results for the NLMS algorithm are given in Figure 
4.5. 

The second and third experiments are devised to compare the convergence per- 
formances of the LMS and NLMS algorithms. The test scenarios for both exper- 
iments are very similar to the first one except that h k = 0.5* with k = 0, 1,2,3, 
[/i y *2 /s] = [0.38'0.56 0.75] for the second one, [/ : f 2 f 3 /a] = [0.38 0.56 0.75 0.57] for 

the third one, and each one has two choices of fi (LMS) for comparison: — 

tr[R vy ] 

and /i 2 = where the total input power trfRyy] is estimated assuming M = L. 

Notice that in the third experiment, two spectrally close sinusoids (/ 2 = 0.56 and 
fi = 0.57) are included in the reference input to create larger eigenvalue disparity 
than that of the second experiment. Results of these two experiments are given in 
Figures 4.6 and 4.7, where we can see that the difference between both algorithms’ 
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(b) time n ; M - 3, L - 4 


(d) time n ; M • 5, L - 4 


Learning curves: lOiog(MSE) of LMS-ANC 



Figure 4.4. LMS-ANC convergence behavior: /x = 0.0093, L = 4; undermodeling 
M = 3 < L, ES = 2.98: (a) Re{w(n)}, (b) Im{w(n)}; overmodeling M = 5 > 
L, ES = 21.60: (c) Re{w(n)}, (d) Im{w(n)}; (e) learning curves: 10 log(J(n)) vs. 
time (ES = 16.75 when M = L). 
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Undermoddling : Re{w(n)} Ov»rmod*Wng: R©{w(n)} 




Und«rmodelling: lm{w(n)} 



Ov®rmod«lling : lm{w(n)} 
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Figure 4.5. NLMS-ANC convergence behavior: fi = 1.0, L = 4; undermodeling 
M = 3 < L, ES = 2.98: (a) Re{w(n)}, (b) Im{w(n)}; overmodeling M = 5 > 
L , ES = 21.60: (c) Re{w(n)}, (d) Im{w(n)}; (e) learning curves: 10 log(J(n)) vs. 
time (ES = 16.75 when M = L ). 


Learning curves: lOlog(MSE) of NLMS-ANC 
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(b) Bme n (•) Bme n 



(c) Bin* n (I) Bme n 


Figure 4.6. Performance comparison for LMS ((a),(b),(c) : p x = 0.0095; (d),(e) ,(f): 
Hi = 0.0119) and NLMS (fi = 1.0), L — 4, reference input SNR = 5 dB, [/x f 2 /a] = 
[0.38 0.56 0.75]; (a),(d): M = 3, ES = 9.67; (b),(e): M = 4, ES = 18.25; (c),(f): 
M = 5, ES = 18.82. 
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undermodelling : M < L undermodeling : M < L 



(a) tima n (d) time n 


correct modelling : M - L correct modelling : M ■ L 



(b) time n (e) time n 


overmodelling : M > L overmodelling : M > L 



(c) time n (f) time n 


Figure 4.7. Performance comparison for LMS ((a),(b),(c): /ij = 0.0073; (d),(e) 
,(f): = 0.0092) and NLMS (ft = 1.0), L — 4, reference input SNR = 5 dB, 

j/i h h A] = [0-38 0.56 0.75 0.57]; (a),(d): M = 3, ES = 12.88; (b),(e): 
M = 4, ES = 29.83; (c),(f): M = 5, ES = 32.97. 
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convergence speeds increases as the eigenvalue spread of increases. This in- 
dicates that the NLMS algorithm’s convergence response is less susceptible to the 
signal’s eigenvalue disparity. 

On the other hand, though the LMS algorithm shows slower convergence through- 
out the simulations, it is able to settle at a smaller value of J( oo), especially in the 
M > L cases. However, increasing p to boost up the LMS algorithm’s convergence 
rate tends to induce larger misadjustment and possible divergence of J(n) sequence 
when the adaptive filter is undermodeling the correlation between both filter inputs. 
Such phenomena can be observed in Figures 4.6(a), 4.7(a), 4.7(e), and 4.7(d). Based 
on these simulation results, with the input signal modeled as a stationary sinusoidal 
process as discussed in Section 4.2, the NLMS algorithm (with p = 1) generally out- 
performs the LMS algorithm in terms of a better tradeoff between convergence speed 
and misadjustment, and stability under a more ill-conditioned signal environment. 

4.4 Adaptive Noise Canceler as Line Enhancer 
As we have addressed before, adaptive noise cancelation requires the presence 
of a reference signal highly correlated with the noise component interfering with the 
signal of interest at the primary input. However, there are several circumstances 
where only one noise-contaminated signal d(n) is available. This problem occurs 
particularly when a broadband signal s(n) is corrupted by a sinusoidal interference 
x(n), and no external reference input free of the signal is available. In such a case, 
the signal d(n) provides its own reference signal y(n), which is taken to be a delayed 
replica of d(n), i.e., y(n) = d(n — A). 

Figure 4.8 shows the block diagram of an adaptive line enhancer (ALE) configured 
via the ANC. Suppose the signal d(n) consists of two signal components: a narrowband 
component x (n) according to (4.12) that has long-range correlations, and a broadband 
component s(n) which tends to have short-range correlations. To design an ALE as 
depicted in Figure 4.8, the delay A is usually selected so that 

E[s(n)s’(n — A:)] ss 0, k > A. 
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reference input 
y(n) = d(n- A) 

= s(n-A) + x(n -A) 

Figure 4.8. Line enhancer via adapti 
input by suppressing the narrowband i 
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Since A is larger than the effective correlation length of the wideband component, 
the delayed replica s(n — A) will be uncorrelated with s(n) in the primary signal, and 
thus the adaptive filter will not be able to respond to this component. On the other 
hand, since the autocorrelation of x(n) does not taper off, the delayed replica x(n — A) 
that appears in the reference input will still be correlated with the narrowband part 
of the primary signal, and the filter will respond to cancel it. 

From our discussion and simulation result in previous sections, the tap weight 
vector in the ALE using the LMS algorithm will approximately converge to the Wiener 
solution when the input signal consists of random-phased sinusoids plus white noise. 
However, the presence of a wideband signal can complicate the statistical character- 
istics of the input autocorrelation matrix R xr . In next chapter, we will investigate 
the performances of LMS-ANC (ALE) in Doppler weather radar applications, where 
the signal environment is very similar to what we propose here. 
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CHAPTER 5 


APPLICATION TO WINDSHEAR DETECTION PROBLEM 

5.1 Weather Detection via Airborne Doppler Radar 
Pulse Doppler radar is used for remote detection of windshears. With appropri- 
ate signal and data processing, this radar can map wind velocity vs. range over an 
antenna scan sector, enabling it to look into storm systems that may contain haz- 
ardous windshears. A form of low altitude windshears known as a microburst has 
been identified as a particular hazard to aircrafts during takeoff and landing [12]. A 
microburst can be qualitatively considered as either wet or dry, dependent on its radar 
reflectivity [52]. Dry microbursts will call for more sophisticated signal processing, 
since the primary sources of reflections, particles of dust and insects, exhibit a much 
lower reflectivity of radar energy. Furthermore, when an aircraft is in low altitude 
flight, such as during takeoff or landing, a portion of the antenna beam is likely to 
illuminate objects on the ground, such as buildings, trees or cars on a freeway. This 
scenario is depicted in Figure 5.1. 

Most of the clutter energy appears around zero Doppler with reference to the 
aircraft ground speed, resulting from strong returns from stationary objects through 
the main beam of the antenna. Additional discrete clutter due to returns from moving 
objects on the ground, or returns from large objects through antenna sidelobes may 
appear at frequencies shifted away from zero Doppler. The large radar cross-section 
of these reflectors on the ground can backscatter enough energy to mask out any 
returns due to weather [53]. This complicates the detection of low-reflectivity weather 
phenomena, where most spectrum-oriented signal processing schemes will operate 
poorly in the presence of a severely low SCR. 




Runway 


Ground 


Figure 5.1. Microburst and wind speed profile. 
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Recently, researchers have suggested that the weather return and the clutter 
return be assumed statistically independent, and setting a signal processing goal of 
decorrelating a nonstationary wideband process (weather) from a nonstationary nar- 
rowband process (clutter) on a per-range cell basis, some preliminary analysis results 
have shown the potential applicability of LMS-based adaptive filters in solving this 
particular signal processing problem [14]. In the next two sections, we will try to fur- 
ther understand the nature of the interfering clutter return that deteriorates weather 
estimation, and to explore weather detection performances of various adaptive and 
passive approaches. 

5.2 Modeling Airborne Doppler Radar Clutter Return 
As elaborated in Chapter 1, the motivation for using adaptive filtering techniques 
in radar clutter rejection is that the separation of two signal processes can be bet- 
ter achieved by exploiting their temporal characteristics. Based on our discussions 
in Chapter 4 relative to the performance of adaptive filters, understanding of the 
interference (or noise) characteristics is more important than knowing the nature of 
the signal. Therefore, before the adaptive filtering techniques can be brought into 
this particular application, the statistical nature of the clutter processes needs some 
investigation. 

Based upon examination of many clutter returns, especially those collected with 
the antenna scan angle kept between ±5 degrees, in terms of their characteristics like 
dominant frequencies, spectral bandwidths, and time-domain fluctuations of average 
power level, we can model the clutter return from a specific range cell as multiple 
sinusoids of constant frequencies, amplitudes, and phases, mixed with complex zero 
mean white Gaussian noise, i.e., 

p 

y( n ) = H A exp(j'27r/,n + <f> { ) + w(n). 

1 = 1 

The signal detection and parameter estimation problems associated with this signal 
model have been addressed in Chapter 2. As the simulation result in Chapter 3 
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testifies, the order-recursive EM frequency estimation algorithm combined with the 
MDL criterion demonstrates excellent performances including high spectral resolu- 
tion, signal detection capability, low SNR threshold, and the statistical characteristics 
possessed by ML estimators. Therefore, this algorithm is used for modeling clutter 
returns. 

Figures 5.2 and 5.3 show two case analyses of typical clutter returns collected 
from the Denver Stapleton airport, with the antenna main beam intercepting the 
ground during a level flight over an interstate highway. Figures 5.4 and 5.5 give similar 
analyses of clutter returns collected from the Philadelphia airport, where the aircraft 
is in a landing approach at about 300 ft above the ground. These spectra show the 
presence of dominant zero Doppler ground clutter and discrete clutter returns due to 
vehicles along the highway. Another two case analyses, as given in Figures 5.6 and 5.7, 
are associated with zero Doppler ground clutter returns collected from the Orlando 
airport. For these analyses, the model order (number of sinusoids) considered for 
the MDL order selection criterion ranges from 1 to 12, and the order p thus decided 
is indicated in each figure. Looking closely at these plots, one can see that our 
algorithm is able to provide a signal model which captures both the spectral mode 
and the temporal characteristics of the radar ground clutter returns. 

5.3 Clutter Rejection via LMS-based Adaptive Filtering 

5.3.1 Experiment for Performance Test 
In order to independently measure the weather detection performances of differ- 
ent signal processing methods, a priori knowledge of both the weather and clutter 
returns should be available. To make this possible, we undertake an experiment as 
depicted in Figure 5.8, where a block diagram shows how real clutter returns can be 
merged with simulated weather returns to test clutter rejection capabilities of adaptive 
filters and fixed-notch filters. The windshear return is regenerated by Research Trian- 
gle Institute’s “Airborne Windshear Doppler Radar Simulation (AWDRS)” program 
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Multiple sinusoids plus complex white noise, ^sequence 



Multiple sinusoids plus complex white noise, Q-sequence 




-20 0 20 
[m/sec] 


Spectrum of estimated clutter signal 



[m/sec] , p - 7 


Figure 5.2. Modeling Denver ground clutter via EM algorithm and MDL criterion, 
p = 7, (NASA test flight-dn4cls5.ml8, frame-77, RC-65). 
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Multiple sinusoids plus complex white noise, l-sequence 



time n, SNR estimate - 24.42 dB 



Spectrum of real clutter return Spectrum of estimated clutter signal 



[m/sec] [rp/sec] , p - 10 

Figure 5.3. Modeling Denver ground clutter via EM algorithm and MDL criterion, 
p= 10, (NASA test flight-dn4cls5.ml8, frame-123, RC-50). 
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Multiple sinusoids plus complex white noise, l-sequence 



Multiple sinusoids plus complex white noise, Q-sequenoe 



Spectrum of real dutter return 



[m/sec] 


Spectrum of estimated clutter signal 



[m/sec] , p - 2 


Figure 5.4. Modeling Philadelphia ground clutter via EM algorithm and MDL crite- 
rion, p = 2, (NASA test flight-f3r3s6.ml4, frame-1036, RC-55). 
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Multiple sinusoids plus complex white noise, Q-sequenoe 



10 20 30 40 50 60 70 80 90 100 

time n 


Spectrum of real clutter return Spectrum of estimated clutter signal 



-20 0 20 -20 0 20 
[m/sec] [m/sec] , p - 7 


Figure 5.5. Modeling Philadelphia ground clutter via EM algorithm and MDL crite- 
rion, p = 7, (NASA test flight-f3r3s6.ini 4, frame-3496, RC-39). 
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Multiple sinusoids plus complex white noise, l-sequence 



Multiple sinusoids plus complex white noise, Q-sequence 



Spectrum of real clutter return 



Spectrum of estimated clutter signal 



[m/sec] , p - 5 


Figure 5.6. Modeling Orlando ground clutter via EM algorithm and MDL criterion, 
p = 7, (NASA test Aight-ol8c43s4.m6, frame-1232, RC-33). 
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Multiple sinusoids plus complex white noise, l-sequenoe 



Multiple sinusoids plus complex white noise, Q-sequence 



Spectrum of real clutter return 


Spectrum of estimated clutter signal 




[m/sec) [m/sec] , p - 7 

Figure 5.7. Modeling Orlando ground clutter via EM algorithm and MDL criterion, 
p = 5, (NASA test flight-o!8c43s4.m6, frame-2394, RC-53). 
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Figure 5.8. Schematic diagram for generation and processing of weather-plus-clutter 
return. 
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[54], assuming the same aircraft operation scenario where the near terminal ground 
clutter return is collected. 

The data submitted for performance test comes from 156 range cells that contain 
clutter returns collected from the Denver Stapleton airport and the correspondingly 
simulated windshear returns, and each range cell has 96 signal samples. To create a 
specific SCR condition, the average power of the simulated weather return is adjusted 
with respect to that of the real clutter return on a range-cell by range-cell basis. 
After this SCR power adjustment, the weather return (signal) and the clutter return 
(interference) are designated as s(n) and x(n), respectively. The weather- plus-clutter 
return, d(n) = s(n ) + x(n), is then submitted for the performance test of different 
adaptive and fixed clutter rejection approaches. 

In our simulation, two adaptive filter configurations, ANC (2nd order) and ALE 
(3rd order, A = 10) implemented with either LMS or NLMS algorithms, are consid- 
ered. We denote them as ANC-LMS, ANC-NLMS, ALE-LMS, and ALE-NLMS. For 
the ANC, the “reference clutter” input y(n ) is generated as a linear combination of 
the clutter return’s delayed replicas, i.e., 

L - 1 

!/( n ) = 51 h k x(n - k) 

k=0 

where denoting H(z) and W 0 (z ) respectively as the z-transforms of hf. and the opti- 
mum Wiener solution w 0 k, h = [h 0 • h.L- i] and the filter’s number of taps M are 

chosen such that H(z ) For the ALE, the number of delays A is chosen so 

as to decorrelate s(n) and s(n — A), i.e., £[s(n)s*(n — A)] ss 0. Remember that the 
ALE is the ANC with its reference input replaced by the delayed primary input, i.e., 

y(n) = d{n — A) = s(n — A) + x(n — A). 

Another clutter rejection filter brought into comparison is a Butterworth second order 
filter with 3 m/sec notch bandwidth centered at zero Doppler velocity. 
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Since the available data record in each range cell contains only 96 samples, these 
LMS-based adaptive filters can not be expected to converge (e.g., refer to Figure 
4.4(e)), though a certain amount of clutter interference can still be suppressed. Fur- 
thermore, experience shows that to reduce the adaptation noise, smaller values of p 
and fi (compared to the upper bounds given in (4.21) and (4.24)) are necessary. In 
order to alleviate this convergence problem and help stabilize the adaptive filters, 
both the ANC and ALE operate on the data twice by retaining the final filter tap 
weights obtained from the first processing as the initial tap weights for the second 
processing. 


5.3.2 Simulation Results 

Figures 5.9 and 5.10 demonstrate the weather detection performance of different 
clutter rejection schemes in terms of their standard deviations of post-processing wind 
velocity estimates, over a wide range of SCR (-25dB ~ lOdB) conditions. In Figure 
5.9, we consider 120 range cells that contain only near zero Doppler ground clutter 
return and weather return, whereas in Figure 5.10, 36 more range cells that have 
both zero Doppler and discrete clutter returns are also taken into account. Based on 
these simulation results, supposing the allowable standard deviation in wind velocity 
estimates is about 2 m/sec, several observations can be noted: 

1. Basically, the presence of discrete clutter does not affect the performance of 
the ANC because the same amount of a priori knowledge regarding discrete 
clutter returns is available, though it significantly degrades that of the fixed 
Butterworth notch filter. However, when the SCR value is above a certain 
level (say, 0 dB), using the Butterworth filter is still an efficient method of 
clutter rejection. 

2. On the other hand, the performance of both ALE-LMS and ALE-NLMS is 
also less sensitive to the discrete clutters than the Butterworth filter’s. This 
can be seen in Figure 5.10 (in comparison with Figure 5.9), where the increase 
in the standard deviation of the ALE windspeed estimate caused by discrete 
clutter is less than 1 m/sec, whereas use of fixed notch filters increases the 
standard deviation by 2 m/sec. 
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3. When discrete clutter is absent, fixed notch Butterworth filters can efficiently 
suppress the zero Doppler ground clutter even with the SCR down to —10 
dB. However, ANC-NLMS, ALE-LMS, and ALE-NLMS can increase the pro- 
cessing gain by at least 7 dB when zero Doppler ground clutter is the only 
hindrance to weather estimation. Considering the effects of the discrete clut- 
ter, adaptive filtering can obtain a processing gain of 7 to 15 dB. 

4. For the ANC, the NLMS algorithm consistently outperforms the LMS algo- 
rithm for the p and /} chosen, whereas ALE-LMS and ALE-NLMS have about 
the same performance. Curiously, for some SCR values, the ALE is able to 
outperform ANC. It may indicate the potential applicability of the adaptive 
line enhancer in windshear detection, when the collection of reference clutter 
returns required of the ANC implementation is technically difficult. 


68 



9 


Standard deviations of wind velocity estimates 



Figure 5.9. Performance comparisons for LMS-based adaptive filters and fixed-notch 
Butterworth filter, considering range cells that contain zero Doppler ground clutter 
and weather returns, (source of clutter: NASA test flight-dn4cls5.ml8). 
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Standard deviations of wind velocity estimates 



SCR, (dB) 


Figure 5.10. Performance comparisons for LMS-based adaptive filters and fixed- 
notch Butterworth filter, considering range cells that contain zero Doppler ground 
clutter, discrete clutter, and weather returns, (source of clutter: NASA test flight- 
dn4cls5.ml8). 
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CHAPTER 6 


CONCLUSIONS 

This work basically deals with three research topics in signal processing and its 
application: (1) signal modeling and parameter estimation, especially for narrow- 
band signals modeled as multiple sinusoids plus white Gaussian noise, (2) LMS-based 
adaptive filter convergence response in this particular signal environment, and (3) 
application of these signal processing techniques in airborne Doppler weather radar. 
Concluding remarks regarding these three areas are presented in the next section. 

6.1 Discussion of Results 

6.1.1 Frequency Estimation and Signal Modeling 

Our focus of investigation is on signals that consist of multiple sinusoids plus 
white Gaussian noise, particularly when only a short record of observations is avail- 
able. Thanks to the simplicity of the complete data model and the monotonically in- 
creasing conditional likelihood function, the EM algorithm is able to solve the problem 
of maximizing the highly nonlinear likelihood function of the signal frequencies. 

Simulation results show that even when the SNR level is very low, our algorithm 
has high frequency resolution capability which most periodogram- or AR-oriented 
approaches fail to possess. The building block of this algorithm is a one-dimensional 
ML frequency estimator, which can be implemented via the Newton-Raphson method 
or FFT, depending on the required numerical accuracy. This indicates the compu- 
tational flexibility of our algorithm. Taking a closer look at Figure 3.1(a), one can 
find another feature of computational parallelism in the sense that all the signal pa- 
rameters (frequencies) can be simultaneously estimated. Further, according to its 
order-recursiveness, for each model order m < p, the initial frequency estimate is 



the optimal (in the ML sense) estimate of previous order (m — 1) plus one additional 
estimate extracted from the model residual z m _j. 

For order selection, we propose the use of Fisher’s T a statistic and the MDL 
criterion tailored for complex data. High probability of correct signal detection of 
these two methods has been observed in Monte Carlo simulation. 

6.1.2 LMS Adaptive Filtering with the Sinusoidal Process 

As an extension of Glover’s work [47] for the case of the complex- valued signal, 
a notch filter implemented by the adaptive noise canceler for suppressing multiple 
sinusoidal interference is investigated through analysis and simulation. Results show 
that a notch filter with deep and narrow spectral nulls indeed can be realized by 
the ANC. Regarding the ANC’s performance response to the sinusoids-plus-noise 
environment, several observations can be made based on our simulation results: 

1. For the LMS algorithm, to maintain stability, the step size p should be small 
compared to the upper bound specified by independence theory, whereas the 
NLMS algorithm with p = 1 generally offers improved performance in terms of 
a better tradeoff between convergence speed and misadjustment, and stability 
in a more ill-conditioned signal environment. 

2. The LMS and NLMS algorithms are capable of approximately tracking the op- 
timum Wiener solution and thus totally rejecting the reference-input-dependent 
component present in the primary input, when the filter’s number of taps is 
greater than or equal to the order of correlation between both inputs. There- 
fore, in most applications where the correlation property associated with the 
signal and noise is unknown, filters of higher order should be employed. 

3. The NLMS algorithm’s convergence response is less sensitive to the eigenvalue 
disparity associated with the signal’s autocorrelation matrix than that of the 
LMS algorithm, though not as significantly as claimed by Slock [28]. It should 
be noted that these nice features associated with the NLMS algorithm only 
come with extra computational cost. 


6.1.3 Windshear Detection Application 
Based upon our analysis and the simulation results, clutter returns have been 
modeled as multiple sinusoids plus Gaussian noise using the EM frequency estimation 
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algorithm and the MDL criterion. The goodness of fit, in terms of the clutter process’s 
temporal and spectral characteristics, demonstrates the robustness of our algorithm. 

Simulation results show that when the SCR level is sufficiently high and only 
zero Doppler clutter is present to bias the estimation of weather information, the 
fixed Butterworth filter is capable of removing the clutter return and thus enhancing 
the weather return. However, when both zero Doppler ground clutter and discrete 
clutter are present in a very low SCR situation, adaptive filtering approaches can 
significantly improve the windshear detection capability. Considering the performance 
of all clutter rejection methods, either adaptive or passive, Figure 5.10 showed that the 
ANC-NLMS stands out as the best approach. Nevertheless, the overall performance 
of the ALE demonstrates a potential applicability in this particular signal processing 
problem, when the implementation of the adaptive noise canceler is impractical. 

6.2 Suggestions for Future Work 

Though our EM frequency estimator produces simulation results that reveal the 
features of the MLE, an analytical result that guarantees the global maximization 
of the likelihood function via the EM algorithm is still not available. Since the EM 
algorithm is potentially capable of solving many MLE problems in signal processing, 
the problems of how to simplify the computation involved for real-time application 
and the analytical verification of its convergence properties in a case-specific fashion 
will be an area for future work. To model clutter returns using the EM frequency 
estimation algorithm and the MDL criterion, our focus is on the narrowband clutter 
data collected with the antenna scan angle kept between ±5 degrees. Clutter collected 
otherwise may have broader bandwidths, and alternative modeling approaches may 
be necessary. 

As for the adaptive filters, though the LMS (NLMS) algorithm is popular due to 
its computational efficiency and ease of implementation, its convergence is slow, and 
the statistical analysis of it for many practical applications, e.g., the realistic signal 
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model proposed in Section 4.2, can be mathematically intractable. The independence 
theory and most analysis results can only provide a picture of this algorithm’s sta- 
tistical behavior from a certain perspective. Furthermore, in a radar application, the 
data available in each range cell is usually quite limited, and the autocorrelation ma- 
trix of the clutter return can be very ill-conditioned. The convergence performance 
of the LMS-based adaptive filters in this type application is an open question. How 
to design adaptive filters with real-time potential, fast convergence, and numerical 
stability, indeed warrants further investigation. 
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APPENDIX A 


ML FREQUENCY ESTIMATOR FOR A SINUSOID PLUS CWGN 


Consider the signal y(n) of a single sinusoid plus complex white Gaussian 
as described in Section 2.2.1, i.e. y{n) = Aexp(j’27r/n) + u;(n), n = 0, 1, • • • , A 

A 

According to (2.2), the MLE of / is given as / = argmax/ jfJ(f) where 
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Given an initial frequency estimate /o obtained from a coarse FFT, the Newton- 



Raphson iterative scheme calculates 
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